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We study theoretically Bose-Einstein condensates with polarized dipolar interactions in anisotropic traps. We 
map the parameter space by varying the trap frequencies and dipolar interaction strengths and find an irregular- 
shaped region of parameter space in which density-oscillating condensate states occur, with maximum density 
away from the trap center These density-oscillating states may be biconcave (red-blood-cell- shaped), or have 
two or four peaks. For all trap frequencies, the condensate becomes unstable to collapse for sufficiently large 
dipole interaction strength. The collapse coincides with the softening of an elementary excitation. When the 
condensate mode is density-oscillating, the character of the softening excitation is related to the structure of the 
condensate. We classify these excitations by linear and angular characteristics. We also find excited solutions 
to the Gross-Pitaevskii equation, which are always unstable. 

PACS numbers: 03.75.Lm, 



I. INTRODUCTION 

Dipolar Bose-Einstein condensates (BECs) have been pro- 
duced using a variety of atoms with magnetic dipoles 0J[4]|. 
In this system the constituent atoms interact via a long-ranged 
and anisotropic dipole-dipole interaction (DDI), which has 
been identified as a source of interesting new effects in the 
degenerate regime, such as supersolidity solitons (9) 

and rotonic excitations fTOlfTTll . 

To realize these novel effects, the DDI must be significant 
compared to the (short-range) s-wave interaction between 
atoms. For 52 Cr (with dipole moment [i m ~ 6yu#) the dipolar 
interaction has been made dominant by using a Feshbach res- 
onance to suppress the s-wave interaction (121 [T3). With the 
recent realization of dipolar BECs of 164 Dy Qi m « lOyUg), 
and 168 Er [4 ] (ji m « 7ju B ) there is now a rich and diverse set of 
systems for experimental investigation. Furthermore, progress 
towards the production of degenerate polar molecules fHifTTl 
with large electric dipole moments promises an exciting future 
in this research field. 

In experiments the dipoles are typically polarized by an ex- 
ternal field, which induces anisotropic (magneto strictive) de- 
formations in the trap 1 18-20] and during ballistic expansion 
l2T1l . Because the DDI has an attractive component (dipoles 
in a head to tail configuration attract each other), an impor- 
tant consideration is in what circumstances the condensate is 
mechanically stable from collapsing to a high density state 
1221 . Experimental USES and theoretical I24H271 studies 
have shown that this stability is highly dependent upon the 
trap geometry. The approach to instability of the conden- 
sate is associated with the softening of elementary excitations 
EHH221- When the condensate is tightly confined in the direc- 
tion that the dipoles are polarized, uniform studies predict that 
high momentum modes will soften fT0ll30lL causing a roton- 
like dip in the dispersion relation (similar to that observed in 
superfluid He EU (32)). In trapped dipolar BECs, the spec- 
trum of elementary excitations is discrete and the excitation 
modes are not momentum eigenstates, however calculations 
have verified that modes with high momentum components 
soften l24ll25ll33U38l. 



Another important feature of dipolar BECs is the emer- 
gence of structured or density-oscillating ground states in 
which the peak density does not occur at the trap center. These 
states have been extensively studied in the cylindrically sym- 
metric case where the condensate can develop a biconcave 
(red-blood-cell) El EH ED or dumbbell (26) shaped density 
profile. Interestingly, it was shown in [ 24 ] that when a bicon- 
cave condensate approaches mechanical instability an angu- 
lar roton (mode with non-zero angular momentum) softens, 
whereas for the normal (non-density-oscillating) condensate 
the mode that softens is purely radial in character (a radial 
roton). 



A study by Dutta et al. 114011 considered the structure of 
the condensate in the more general case of a fully anisotropic 
trap, and found an array of different density-oscillating ground 
states. That study motivates the work we report here, where 
our aim is to investigate the nature of the excitations in the 
fully anisotropic trap, particularly those that soften in the ap- 
proach to mechanical instability. Indeed, our work extends 
that of l24l to show that a variety of rotonic modes, which 
can be broadly classified by angular or linear characteristics, 
emerge in dipolar BECs. These rotons usually reveal some 
character of the underlying condensate, particularly if it is in a 
density-oscillating state. We note that the nature of the soften- 
ing modes has been suggested as a probe of collapse in l25l . 
During our study we also found that our ground state phase 
diagram was inconsistent with the results in Ref. l40l . 



The remainder of the paper is organised as follows: in Sec- 
tion|ll| we outline our model and methods of solution; in Sec- 
III ApII B" we discuss the structure and stability of dipo- 



tions 



lar BEC ground states; in Section [nrC|we c ompare our results 
to those in Ref. l40ll ; and in Section |III D| we discuss the na- 
ture of the elementary excitations whose softening accompa- 



nies the collapse. We conclude in Section IV 
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II. FORMALISM 

We consider a BEC of N bosons trapped in an external po- 
tential 



Vext(r) = (qj 2 x x 2 + oj 2 y 2 + co 2 z z 2 ) , 



(1) 



where the trap frequencies 0J x , y , z ma Y be distinct. We take the 
atoms to have magnetic dipole moment of magnitude \i m po- 
larized by an external field along the z direction. The resulting 
DDI potential between atoms is of the form PHI 



V dd (r) = 



jdojij 1 - 3 cos 2 
47r r 3 



(2) 



where //o is the permeability of free space and is the angle 
between r and and the z axis. In general the effective low 
energy interaction for dipolar atoms also includes a contact 
interaction, i.e., 



47777 fp- 

Vint(r) = — 6(r) + V dd (r), 



111 



(3) 



where a s is the s-wave scattering length. However, here we fo- 
cus on the purely dipolar case with a s = 0, as can be arranged 
in experiments by use of Feshbach resonances (e.g. see (T3)). 



(quasi-particle) excitations of the condensate, described by the 
Bogoliubov-de Gennes (BdG) equations 



where 
L 



The operators xi and xi 



Xif(r) =N f JrV*(r0y in t(r-r0/(r0^(r), 
Xifix) =Nf Jr^CrOVintCr-rO/CrO^r), 



(8) 



(9) 



(10) 
(11) 



describe the exchange interactions between the condensate 
and excited modes, and Q is the projection onto the subspace 
orthogonal to the condensate: 

G/(r) = /(r) - <A(r) J drV*(rW). (12) 

The method of solution of the BdG equations ^ is explained 
in Appendix|A2|. 



A. Condensate mode 

The equilibrium unit-normalized mode i/f(r) for the conden- 
sate is determined by the non-local Gross-Pitaevskii equation 
(GPE) (13: 



l*Kr) = [H +N J Jr , y int (r-r , )|<A(r , )| 2 ]<A(r), 
= Xcp<A(r), 

where ji is the condensate chemical potential, 

h 2 



2m 



V z + VextW, 



(4) 
(5) 

(6) 



is the single particle Hamiltonian, and the last (convolu- 
tion) term in Eq. ^ describes the direct (Hartree) interac- 
tion of condensate atoms with themselves. We solve for self- 
consistent solutions of Eq. ( [j) us ing a Newton Krylov algo- 
rithm, described in Appendix A 1 



B. Elementary excitations: Bogoliubov-de Gennes equations 

Linearizing about the GPE solution if/(r) the time- 
dependence of the condensate l42l is given by 

¥(r, t) = {<A(r) + -±= ^ [cjUjWe-™* + c*v*(r)^] }e^\ 

(7) 

where the c/s are constants (e.g. see 1431 ). The modes 
{uj,Vj}, with respective frequencies ojj, are the elementary 



III. RESULTS 

Here we consider a harmonic trap which is, in general, fully 
anisotropic, i.e., the three trap frequencies are different. In 
our results we use cj y and l y = ^Jh/maj y as the reference fre- 
quency and length scale, and characterize the trap by the two 
anisotropy parameters 



K - OJy/OJy, 



(13) 
(14) 



It is important to note that even if A x = A z = 1 (spherical trap), 
the system is only cylindrically symmetric about the z axis 
because the dipoles are polarized along this direction. Thus, 
whenever A x ^ 1 the system breaks cylindrical symmetry, and 
requires the fully three-dimensional numerical solution. We 
also introduce the dimensionless interaction parameter 



D = 



4tt h 2 l' 



(15) 



which is the same as that used in (24l [29) in the cylindrically 
symmetric limit. Any ground state is uniquely specified by 
the parameter tuple (A x , A Z ,D). 

It is worth noting that the results presented in [40 ] made 
the unconventional choice of dipoles polarized along the y- 
direction and used the x-oscillator parameters to define their 
units. Our results can be compared to theirs by cyclically per- 
muting the coordinates, i.e. A x m <-> A 



their 



their 
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FIG. 1 . (Color online) (a) Stability plot of a dipolar condensate in an 
anisotropic trap as a function of the anisotropy parameters A x and A z 
and dipole strength D. In the dark- shaded (magenta) region there ex- 
ist stable density-oscillating condensates (with the maximum density 
away from the center of the condensate). The light- shaded (cyan) re- 
gions are normal ground states with maximum density at trap center. 
The slice in the A z direction, corresponding to cylindrical symmetry 
A x = 1 , is displaced from the axes to improve visibility and repro- 
duces Fig. 1 of Ref. (24). (b)-(e) Stability plots for constant A z =6 
(b), 6.5 (c), 7 (d) and 8 (e) with the type of state labelled I-V accord- 
ing to the classification given in Sec. |IIIB| (f)-(g) Stability plots for 
constant D = 15 (f) and 30 (g). 



A. Ground state stability 

Our key results summarizing the stability diagram are 
present in Fig. [TJa). The shaded regions in this plot indicate 
where a stable ground state can be located, with dark shad- 
ing used to indicate the region in which a density-oscillating 



state occurs (i.e., where the condensate does not have its peak 
density at trap center). Our results for the case of cylindri- 
cally symmetric confinement A x - 1 correspond to those in 
Ref. 1 24], including that the density-oscillating states are bi- 
concave. Dipolar stability is highly dependent upon geome- 
try, as is revealed in Fig.[TJa). A general trend of Fig.[TJa) is 
that stability increases as A z increases and A x decreases. This 
can be understood because the DDI is attractive for dipoles 
in a head-to-tail configuration (z separation) and repulsive 
for dipoles in a side-by-side configuration (xy-plane separa- 
tion). Thus increasing A z (tightening z confinement) reduces 
the number of dipoles in the destabilizing attractive configura- 
tion, while decreasing A x (loosening x confinement) increases 
the number of dipoles in the stabilizing repulsive configura- 
tion. 



B. Ground state types 




FIG. 2. Different types of ground state density profiles. Density slice 
of the condensate in the xy-plane, focusing on central region where 
density oscillations can occur. Solutions give examples of the var- 
ious condensate types (I-V) and are obtained for the parameters I 
(A x , A z , D) = (0.65, 9, 89.4); 11(a) (0.85, 6, 20.3); 11(b) (0.9, 6.5, 26.3); 
III ( 0.6,6 .5,50.1); IV (0.65,6,37.8); V (1,6.5,24). See text in 
Sec. |IIIB| for a discussion of the ground state types. 

We find that a density-oscillating condensate can exhibit a 
range of different shapes, however in the region of interest the 
Z confinement is sufficiently tight that the non-trivial density 
features occur in the xy-plane. Following [40] we categorize 
the ground state solutions by the labels I-V as follows, with 
reference to examples in Fig.[2j 
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type-I Normal (non-density-oscillating) condensate with peak 
density at the origin. 

type-II Density-oscillating condensate with two-peaks in the x- 
direction. This can either be a simple two-peaked struc- 
ture [e.g. Fig. [2111(a)], or two peaks with a biconcave 
crater [e.g. Fig7|2] 11(b)]. 

type-Ill Density-oscillating condensate with two-peaks in the y- 
direction. Like type-II, this case could be a simple two- 
peaked case (not shown) or with a biconcave crater. 

type-IV Density-oscillating condensate with four peaks (two 
peaks along both the x and v-directions). 

type-V Density-oscillating condensate that is biconcave (with 
no additional peaks). 

In Fig. [3jg) we present a projection (bird's-eye view) of 
the stability diagram [Fig. [TJa)] onto the A X A Z -plane, indicat- 
ing the different types of ground states which occur. Nor- 
mal (type-I) condensates are found in the light shaded region 
for any value of D. Dark shaded regions indicate if a stable 
density-oscillating state exists for any value of D in that trap 
geometry. The darkly shaded region is sub-divided accord- 
ing to the type of density-oscillating state (types II- V) at the 
largest dipole strength D for which a stable ground state exists 
at that value of A x , A z . 

Contours of a generalized anisotropy A = oj z / ^oj x oj y are 
shown in Fig. |3jg). This parameter is a measure of the con- 
finement along the direction that dipoles are polarized rela- 
tive to the geometric mean confinement in the xy-plane (mo- 
tivated by the discussion of trap effects on stability at the end 



of Sec. Ill A). We note that A tends to qualitatively charac- 
terize how some boundaries of the density-oscillating region 
develop as the trap geometry changes. 

Figures [T] and [3jg) reveal the dominant role that the type-II 
and III states have in the density-oscillating region. Gener- 
ally type-II states are favored at lower values of A z and D and 
type-Ill states at higher values of A z and D. Intricate struc- 
tures (lobes) in the shape of the density-oscillating conden- 
sate region arises in parameter regimes where both conden- 
sate types are present [e.g. Figs.[TJb) and (c)]. Type-IV and V 
condensates are less common, with type-IV states emerging in 
the transition between the type-II and III regions, and type-V 
states occurring near cylindrical symmetry (i.e. A x « 1). As 
the system becomes more anisotropic [A x < 0.6, see Fig.[3jg)] 
the type-II and III sub-regions become distinct in A x A z -space 
(i.e. separated by a normal type-I region) and appear to emerge 
as two separate branches. We note that our results for A x < 1 
can be mapped onto solutions for A x > 1 by exchanging x 
and y coordinates (e.g. so that type-II states become type-Ill 
etc.) and scaling the interaction parameter. This reveals that 
the type-II and III branches cross at A x = 1 . 



C. Comparison to Dutta et al. (40) 

The parameter regime in Figs.[lland|3]is similar to that stud- 
ied by Dutta et al. [40 ]. In Fig.|4[a) we present the data cor- 
responding to Fig. 2 of 1 40 ]. Our results disagree with theirs 




FIG. 3. (Color online) (a)-(f) Condensate density slices in the 
xy-plane for the parameters: (a) (A X ,A Z ,D) = (0.8,6.5,35.1); (b) 
(0.9,6.5,26.3); (c) (1,6.5,24); (d) (0.6, 6.5,50.1); (e) (0.65,6,37.8); 
(f) (0.85,6,20.3); (g) Projection of upper surface (large D) ground 
state types in Fig. [TJa) onto the /l^-plane. The density-oscillating 
(dark- shaded) region is divided according to the type of ground state 
therein. The light-shaded (cyan) regions have no stable density- 
oscillating ground states for any value of D. For reference, contours 
of constant A-cx) z l ^io x co y are shown (thick lines). 



in the following significant ways [44]: 1. The results in Fig. 2 
[ 40 1 predict stability to higher dipole strengths than our results 
at given trap aspect ratio. 2. We do not find density-oscillating 
states for the same parameters they report. For the density- 
oscillating region in Fig.|4ja), we only find two-peaked states 
of the kind shown in Fig. [4jt>), and not four-peaked states 
found in l40l . We also note that their density-oscillating re- 
gion extends over a broader range of A x values than ours. 3. 
We always find our density-oscillating states to be even with 
respect to x-, y-, and z-reflections, whereas in Fig. 4 of l40l 
ground states are reported which break this symmetry. 

The origin of these differences is not clear to us. The imagi- 
nary time algorithm used in HOI to locate ground states is only 
briefly discussed, although they reported that the results were 
dependent on the initial (random) states used. Our algorithm 
(as outlined in Appendix | A 1| ) is built around robust optimiza- 
tion techniques and has been carefully checked against other 
approaches (e.g. the cylindrically symmetric results reported 
in (29)). We also employ a spherical cutoff interaction poten- 
tial, which has been shown to improve accuracy on finite grid 
calculations by minimizing aliasing effects of the long-ranged 
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FIG. 4. (a) Stability plot of a dipolar condensate as a function of 
A x , with /i z = 5.5. Two-peaked ground states are indicated by the 
dark (magenta) region. The dotted line indicates the plotting range 
of the energy of unstable four-peaked solutions shown in (d). (b) 
Stable two-peaked solution for (A X ,A Z ,D) = (0.65,5.5,20). (c) Un- 
stable four-peaked GPE excitation for (0.85, 5.5, 20). (d) Energy per 
particle of type-I ground state (solid line) and unstable, four-peaked 
(type-IV), GPE excitation (dotted line) as a function of D for the 
same trap parameters as in (c). 



interaction l29l . 

Another important feature of our work is that we confirm 
the stability of our solutions by performing a BdG analysis 
of the excitations. For a stable ground state all the quasi- 
particles have real energies. However, it is possible to find sta- 
tionary solutions of the GPE for which some excitations have 
an imaginary frequency. In this case these excitations would 
grow exponentially in time and the ground state is dynami- 
cally unstable. The importance of this is revealed in Fig.[4jc), 
where we show a 4-peaked solution of the GPE equation that 
we obtained at an interaction strength well-above the stabil- 
ity boundary. This solution has excitations with imaginary 
eigenvalues, so is dynamically unstable. In Fig.Qd) we com- 
pare the energy per particle [Eq. ( A3 >] of the 4-peaked solu- 
tion against that of the dynamically stable 2-peaked solution 
(which ends at the stability boundary, D « 10). Interestingly, 
the 4-peaked solution exists for dipole strengths where the 2- 
peaked solution exists, but is of much higher energy. These 
results demonstrate that BdG solutions are vital in determin- 
ing stable condensate states. 




FIG. 5. (Color online) Bogoliubov excitation spectrum as a func- 
tion of dipole strength Z), for A x = 0.85, A z = 6. The softening 
mode is labelled by the blue dashed line. The colors represent the 
parity along the x-direciton (light/red=even), (dark/black=odd) . For 
these parameters the condensate is a type-I state up until near col- 
lapse where it changes to a type-II state [e.g. see the relevant phase 
diagram Fig.[TJb)] 



observe that a mode, which has a relatively high excitation 
frequency at low values of D, decreases its frequency as D 
increases and eventually approaches zero at the point of insta- 
bility. Because this mode tends to have short wavelength fea- 
tures it is referred to as a rotonic excitation. For definiteness 
in this section we will refer to the lowest energy such mode 
for D values close to instability as the roton mode and will 
label it as the j = 1 quasi-particle [i.e. mode {u\(r), vi(r)}]. 

We will study the properties of these rotons in the same gen- 
eral parameter regime that we used to study the ground states 
i.e. with oj 7 



D. Bogoliubov excitation spectrum: roton modes 



in Sec. Ill A i.e. with oj z » oj x ,oj y . In this regime the roton 
mode is structureless in the z-direction and exhibits structure 
in the xy-plane. Indeed, the analysis of a uniform quasi-two- 
dimensional dipolar condensate confined tightly along the z 
direction predicts that the roton modes lie in the xy-plane with 
a characteristic momentum set by the z-confinement length 
scale fTTll . The study of roton modes in a cylindrically sym- 
metric pancake trap (oj z » oj x = oj y ) showed that its structure 
revealed properties of the condensate state l24l . Noting that 
for the cylindrically-symmetric trap quasi-particles are eigen- 
states of the angular momentum operator L z with eigenvalue 
m z , two cases of roton modes were found: when the con- 
densate was in a normal (type I) ground state a radial roton 
emerged with m z = 0. When the condensate was in a bicon- 
cave (type V) state an angular roton with \m z \ > emerged. 

Our primary concern here is to examine the structure of the 
roton modes in the fully anisotropic trap, and the relationship 
this has to the various types of condensate ground state. We 
begin by introducing the techniques we use to visualize and 
characterize the roton excitations. 



In this section we examine the properties of the quasi- 
particle excitations of a dipolar condensate. Our particular 
interest is in the modes that soften and cause the condensate to 
become dynamically unstable as the dipole strength increases. 
For example, in Fig. [5] we show the excitation spectrum of a 
dipolar condensate as a function of the dipole strength. We 



1. Density fluctuation 

It is useful to consider the effect that the roton, when excited 
to have a small coherent amplitude A = \A\e 1 ^ with |A| <?c 1, 
has on the condensate. In this case the dynamics of the total 
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FIG. 6. (Color online) Bogoliubov excitation spectrum as a func- 
tion of dipole strength D for (a) A x = 0.65, A z = 9; (b) A x = 0.85, 
A z = 6; (c) A x = 0.9, A z = 6.5; Spectra color-coded to indicate (p 2 x ). 
Contours of roton density fluctuations (6n\) in the xy-plane for (Rl) 
D = 20.3 [with type-I condensate]; (R2) D = 26.3 [type-II]; (R3) 
D = 89.4 [type-II]; and trap parameters corresponding to (a), (b), 
and (c), respectively. Thick lines indicate the nodal lines. 



density [i.e. setting c\ — > A and all other Cj set to zero in 
Eq. ([7])] is given by l45l 

n(r, t) = N\ifs(r)\ 2 + 2\A\ cos(a^ - cp) Vtf^(r) Sm(r), (16) 

to linear order in u\ and V\. We have introduced 

Snjir) = Uj(r) + v 7 (r), (17) 

as the density fluctuation associated with the j-th quasi- 
particle. To visualize the roton mode we plot contours of 
6n\ (x) (see Figs. [6] and [7]). 

2. Roton characterization 

We also consider how to generalize the qualitative descrip- 
tion of the roton modes (e.g. radial and angular rotons of l24l ) 
to the anisotropic trap. Because the excitations are not (in 
general) eigenstates of L z this characterization cannot be per- 
formed by inspection of the m z value of the relevant mode. 
Thus, we propose to characterize quantitatively the rotonic 




FIG. 7. (Color online) Bogoliubov excitation spectrum as a function 
of dipole strength D for (a) A x = 0.6, A z = 6.5; (b) A x = 0.65, A z = 6; 
(c) A x = 0.8, A z = 6.5; Contours of roton density fluctuations (dni) 
in the xy-plane for (R4) D = 50.1 [with type-Ill condensate]; (R5) 
D = 37.8 [type-IV]; (R6) D = 35.1 [type-IV]; and trap parameters 
corresponding to (a), (b), and (c), respectively, (also see Fig.[6j. 



modes by computing the expectations of various operators, as 
follows: 

f Jrw*(r)Ow/(r)+ f dr v*(r)Ov ; (r) 
(0)j = ~ l r —r f J —i ' ( 18 ) 

with O being p 2 x , p 2 , or L 2 (where p x - -ihj^. etc.). The 
results of this analysis for the roton modes are given in Ta- 
ble |lj and the {p 2 ) values for all quasi-particles are used to 
color-code spectra in Figs. [6] and [7] We note that Wilson et 
al. lf33l used a similar procedure to assign a momentum to 
each quasiparticle according to pj - ^j{p 2 x )j, and allow them 
to approximately extract a dispersion relation in the cylindri- 
cally trapped gas. We interpret these expectations as follows: 

• Angular character. We take (L 2 ) > 4h 2 to define the 
angular characteristic of a roton mode. Note: We take 
an (L 2 ) value consistent with \m z \ = 2 to define the an- 
gular characteristic because this is the minimum value 
of angular momentum found for an angular roton in the 
cylindrically symmetric case. 
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TABLE I. Roton mode properties for various types of ground states. The small ticks for bi-concavity and angular character qualities of R4 
indicate that these are marginal. Excitations for the two cylindrically symmetric cases R7 and R8 are not shown in figures here, however R8 
corresponds to Fig. 2(Ic) of |24|. Note that there are two degenerate roton modes for R7 with \m z \ > 0, and both have identical squared linear 
momentum and angular momentum expectations. 



• Linear character. When {p 2 ) > (p 2 ) the roton mode os- 
cillates more rapidly along the x direction and we refer 
to the mode as a linear roton along x. Similarly, when 
(Py) > (Px) we nave a Amear roton along y. A radial ro- 
ton has no directional preference, i.e. has (p 2 ) « {p 2 x ). 

Note that these characteristics are not exclusive, for example, 
the analysis R6 in Table [I] reveals a roton that has both linear 
and angular characteristics. We emphasize that our character- 
istics agree with the terminology adopted by Ronen et al. l24l 
for the cylindrically- symmetric case (e.g. cases R7 and R8 in 
Table [I]): R7 has angular character (termed an angular roton in 
l24l ) for a (type-V) biconcave condensate. R8 is a radial ro- 
ton for a cylindrically symmetric normal (type-I) condensate. 
We note that in the cylindrically symmetric case we must have 
(p 2 ) = (py), thus there can be no linear character. 



3. Roton analysis 

In Fig. [6] we give some examples of roton modes for normal 
(type-I) condensates and density-oscillating condensates with 
peaks along the x-direction (type-II). The quantitative analy- 
sis and characterization of these is given as results R1-R3 in 
Table |T| For these examples the roton varies rapidly along the 
x-direction. More generally, the spectra shown in Fig. [6] re- 
veal that, of the low-energy quasi-particle modes considered, 
the modes that most rapidly descend as D increases are those 



with the largest values of {p 2 ). We also note that the roton 
for the two type-II states considered differ in their effect (as 
a density perturbation) on the condensate: in one case [Fig. [6] 
R2] the roton causes the two peaks of the condensate to os- 
cillate out-of-phase, while the other case [Fig. [6]R3] causes 
these peaks to oscillate in-phase (also see l45l ). In contrast, 
for the normal ground state the density fluctuation of the roton 
is strongest at trap center [Fig.[6]Rl] where the ground state 
has its peak density. 

In Fig. [7] we consider cases in which the condensate has 
structure in the j-direction, i.e. type-Ill and type-IV conden- 
sates. The quantitative analysis and characterization of these 
is given as results R4-R6 in Table [I] Interestingly, the emer- 
gence of y peaks does not mean that the roton will be most 
rapidly varying along the j-direction, e.g. as revealed in the 
momentum expectations of R4 in Table [l| However, we do 
observe that the additional y structure tends to emerge in the 
roton (see 6n\ plots in Fig. [7]), and is accompanied by an in- 
crease in (p 2 ). The spectra for these states in Fig. [7] indicate 
that modes with the highest values of {p 2 ) do not preferen- 
tially go soft first. The density fluctuation maxima associated 
with these rotons generally occur at locations corresponding 
to the peak density of the condensate, and dynamically causes 
the peaks oscillate in various ways (e.g. see B31 ). 

A general trend we see is across all states considered is that 
when the condensate has some biconcavity in its structure (as 
discussed Sec. |IIIB }, that the roton tends to exhibit angular 
character (see Table l|. 



IV. CONCLUSIONS 



In conclusion, we have mapped the stability and structure 
of anisotropically trapped dipolar BECs in parameter space, 
and explored the occurrence of density-oscillating conden- 



sates. The parameter regime we have examined is dominated 
by ground states with two peaks along the x or y direction. 
We observe these ground states to form clear sub-regions of 
parameter space, which branch off in the highly anisotropic 
case (A x <*c 1). 

Collapse instability is associated with the softening of a ro- 
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ton excitation, which in anisotropic traps we typically find to 
have a linear character with it being most highly excited in the 
direction of weakest trap frequency. The softening roton mode 
has characteristics associated with the ground state structure, 
e.g., if the ground state has biconcave character then the ro- 
ton tends to have an angular character; we also find that the 
density fluctuations of the roton mode coincide with the peaks 
of the condensate density. By increasing the dipole strength 
slowly until the system is dynamically unstable we would ex- 
pect that this roton structure will be revealed in the collapse 
dynamics (e.g. see l25l ). 

The linear character of the roton mode suggests that the crit- 
ical velocity for breakdown of superfluidity will be anisotropic 
in the xy-plane. This differs from the anisotropic superfluidity 
predicted for pancake trapped dipolar condensates in 1371 , in 
which the anisotropic arises from tilting the dipole polariza- 
tion axis into the xy-plane. It will be interesting to explore 
the nature of the rotons in a fully anisotropic trap under much 
tighter z-confinement. 

As most experiments in dipolar condensates (3] 01 \13\ [22l 
[231 l46l are conducted in anisotropic traps, we believe our 
more general study of condensate and roton structure will as- 
sist in experiments planning to measure roton properties. 



ACKNOWLEDGMENTS 

This work was supported by the Marsden Fund of New 
Zealand contract UOO0924. 



Appendix A: Brief outline of numerical methods 

In this appendix we briefly discuss the approach we use for 
solving the GPE and BdG equations. Our approach is similar 
to that outline in Ref. l29l. We will work in the dimensionless 



units that were introduced in Sec. Ill In these units the single- 
particle operator and DDI potential take the form 



3 cos 2 8 



V D (r) = D- 

respectively, where we have set NV[ nt 
elude the condensate number, N. 



(Al) 
(A2) 

Vd to explicitly in- 



and n u (r) = \i// u (r)\ 2 . To avoid having to deal with the the con- 
densate's normalization constraint we have allowed the con- 
densate orbital to be unnormalized (denoted if/ u ) and explicitly 
scaled out the normalization dependence in E [if/ u ] using the 
normalization functional N = f \if/ u (r)\ 2 dr |471 . This means 

that the un-normalized (i// u ) and normalized {if/ = i// u / ^fN) 
orbitals give the same energy in Eq. dA3j. We also note 



that ®£>(r) can be efficiently evaluated using the convolu- 
tion theorem D (r) = ^F -1 [V/)(k)/T M (k)], where T is a three- 
dimensional Fourier transform, h u (k) = T{n u (Y)} and Vz)(r) = 
T{V D (X)} (which can be evaluated analytically EU (29)). In 
our code T is implemented using forward and reverse fast 
Fourier transform (FFT) algorithms. The FFT implicitly treats 
the system as a 3D lattice of condensates of period 2R, where 
our numerical grid is cubic with extent [-R, R] in the x, y, z 
directions and our results here use R = 8 with 64 points in 
each direction. We adopt a radial truncation of Vz)(r), such 
that Vz)(r) = for r > R, which prevents its overlap with the 
condensate 'copies', giving improved numerical accuracy of 
®£>(r). The Fourier transform of this truncated potential has 
the analytic form |29l 



1 | 3 cos(Rk) 3 sin(Rk) 



R 2 k 2 R 3 k 3 
where a is the angle between k and the & z -axis. 



(3 cos 2 a - l), 

(A5) 



Our general procedure discussed here is similar to that pre- 
sented in [29]. However, we differ from that most signifi- 
cantly in that instead of using a conjugate gradient technique 
we use a Newton-Krylov algorithm l48l to solve for the con- 
densate. This algorithm iteratively finds zeros of the residual: 
dE/dif/* u = [£gt>i/s u ~ ^u] I VVV, where 



£ GP = H + N- l ® D (v), 



(A6) 



and// = f dri//* u £ G ?\J/ u /N. 



2. Solution of Bogoliubov-de Gennes equations 



1. Newton Krylov method of obtaining stationary states 

GPE solutions can be found by minimizing the energy func- 
tional: 



E[l//u] = h f dx ^ u 



IN 



i/su, (A3) 



where 



(D D (r) = J dv' V D (r - r'K(r'), (A4) 



We formulate the BdG matrix £ [Eq. d9l)] in the GPE ba- 



sis, i.e., eigenstates of Eq. (A6), which we calculate directly 
on our numerical grid using the Arnoldi algorithm provided 
by the MATLAB routine 'eigs'. The integrals involving the 
dipolar interaction potential are done analogously to those de- 



scribed in Sec. II A We diagonalize the resulting BdG ma- 
trix using the MATLAB routine 'eig'. For most cases studied 
in this paper, we find 400 basis vectors sufficient for conver- 
gence of the excitation energies cou however, bases of more 
than 1200 vectors are required for condensates tightly trapped 
in the z direction (e.g., for A z = 9). 
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